In [1]:
import pandas as pd
import numpy as np
from pandas import read_csv
from matplotlib import pyplot
import matplotlib.pyplot as plt

Задание 9:

Для заданной страны необходимо построить прогнозы на ближайшие 5 лет показателей деловой активности (индекс промышленного производства, индекс деловой активности, ВВП в постоянных ценах, уровень безработицы). Использовать месячные данные. В качестве возможных факторов использовать: фондовый индекс, уровень процентных ставок, отраслевые индексы цен, курс национальной валюты к доллару США, уровень цен на основные сырьевые товары (нефть и золото), дефлятор ВВП, индекс потребительских цен и основные денежные агрегаты (денежная база, денежная масса (M1, M2)).

Исследовать, улучшают ли выбранные факторы прогноз показателей деловой активности в сравнении с одномерной моделью временного ряда?

Для нашего исследования деловой активности мы решили использовать следующие данные:

Какие данные используются?

  • Индекс деловой активности в производственном секторе (PMI) США от ISM
  • Композитный индекс деловой активности (PMI) в США от Markit
  • Фьючерс на нефть Brent
  • Фьючерс на золото
  • Уровень безработицы
  • U.S. Gross Domestic Product
  • S&P 500 (SPX)
  • Объем производства %
  • Курс Доллара к Евро

Среди показателей деловой активности мы выбрали два похожих по содержанию показаткелей - индексы деловой активности в производственном секторе и композитный индекс деловой активности.

Расскажем немного о них:

PMI — это макроэкономический показатель, позволяющий судить о стабильности экономики. Во всем мире для оценки состояния экономики используют индексы деловой активности. Наибольшее распространение они получили в США, Китае и странах Еврозоны. Эти показатели имеют большое влияние на экономические процессы, но главным образом на рынок ценных бумаг.

Индекс PMI рассчитывается на основании опроса топ-менеджеров ключевых производственных предприятий. Ответы даются на основании сравнения существующих показателей со статистикой истекшего периода (календарный месяц). Значение PMI, которое представляет собой число от 0 до 100, определяет собой оценку изменения Валютного курса страны, Темпы роста ВВП, стоимость ценных бумаг и других ведущих показателей.

Разновидности индекса PMI

  • Активность в производственном секторе
  • Активность в сфере услуг
  • Активность в непроизводственном секторе
  • Композитный индекс

В своем проекте мы решили взять для анализа и предсказания два наиболее популярных индекса- композитный и произвеодственный в качестве предсказываемых переменных.Именно эти индексы наиболее чуствительны к экономическим колебаниям и изменениям в экономике стран. Мы также выбрали в качестве объясныющих переменных основыной индекс S&P 500, который включает в себя акции наиболее крупных Американских компаний, большинство менеджеров которых и определяют индекс PMI. Мы также собрали данные по ВВП страны и основные котировки на сырьевом рынке и рынке драг.металлов: Brent и Золото, а также ИПЦ США.

Цель исследования :

В кругах финансовых аналитиков ходят разные мнения на счет направления влияния индекса деловой активности и определяющих факторов этого индекса.Кто-то говорит, что оценка менеджерами "перспектив" влияет на ВВП, акции и иные финансовые индексы, кто-то считает наоборот.В данном проекте мы попытаемся проанализировать временные ряды, определить их свойства и построить прогноз на ближайшие 5 лет индексов PMI.

Загрузка и визуализация данных:

In [2]:
import matplotlib.pyplot as plt
import plotly.express as px
from datetime import datetime
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
df['month_year'] = pd.to_datetime(df['Date']).dt.to_period('M')
#df=df.sort_values(by='Date')
fig = px.line(df, x='Date', y='PMIcomp')
fig.add_hline(y=50.0,line_dash="dash", line_color="red")
fig.show()
fig = px.line(df, x='Date', y='Vol')
fig.show()
fig = px.line(df, x='Date', y='PMItot')
fig.add_hline(y=50.0,line_dash="dash", line_color="red")
fig.show()
fig = px.line(df, x='Date', y='SP')
fig.show()
fig = px.line(df, x='Date', y='GDP')
fig.show()
fig = px.line(df, x='Date', y='Gold')
fig.show()
fig = px.line(df, x='Date', y='Brent')
fig.show()
fig = px.line(df, x='Date', y='Unemployment')
fig.show()
fig = px.line(df, x='Date', y='USD-EUR')
fig.show()
#df.PMIcomp = df.PMIcomp.astype(int)
#df.PMItot= df.PMItot.astype(int)
# Everything in one
fig = px.line(df, x='Date', y=['Vol','USD-EUR','SP','GDP','Brent','Gold','PMItot','PMIcomp'])
fig.show()

Анализ на стационарность

In [3]:
from statsmodels.tsa.arima.model import ARIMA 
import statsmodels.tsa.stattools 
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
alpha=0.05
#perform augmented Dickey-Fuller test

result = adfuller(df['SP'], autolag='AIC')
print(f'ADF Statistic: {result[0]}')
print(f'p-value: {result[1]}')
for key, value in result[4].items():
    print('Critial Values:')
    print(f'   {key}, {value}') 
    
if result[1] < 0.05:
    print('Ряд стационарен')
else:
    print('Ряд не стационарен')
ADF Statistic: -2.0456393951002196
p-value: 0.2669270483950142
Critial Values:
   1%, -3.4626576734812318
Critial Values:
   5%, -2.8757444215841326
Critial Values:
   10%, -2.5743412314098753
Ряд не стационарен
In [4]:
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['GDP'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-6.426876181767085 1.735758020949155e-08 Ряд стационарен
In [5]:
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Gold'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-2.837032739325091 0.0531755131923 Ряд не стационарен
In [6]:
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Brent'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-1.2594868454127226 0.6474974371691998 Ряд не стационарен
In [7]:
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['PMItot'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-6.138075978190617 8.095526440642527e-08 Ряд стационарен
In [8]:
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['PMIcomp'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-4.033271313743218 0.001245971056526383 Ряд стационарен
In [9]:
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['USD-EUR'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-1.3469781779092913 0.6074153069855851 Ряд не стационарен

KPSS -test. Так как индексы деловой активности строго ограницены [0-100] и колеблются около 50 - мы должны провести иной тест на стационарность, диагностирующий наличие либо отсутствие тренда.

In [10]:
from statsmodels.tsa.stattools import kpss
import warnings
warnings.filterwarnings("ignore")
def kpss_test(series, **kw):    
    statistic, p_value, n_lags, critical_values = kpss(series, **kw)
    # Format Output
    print(f'KPSS Statistic: {statistic}')
    print(f'p-value: {p_value}')
    print(f'num lags: {n_lags}')
    print('Critial Values:')
    for key, value in critical_values.items():
        print(f'   {key} : {value}')
    print(f'Result: The series is {"not " if p_value < 0.05 else ""}stationary')

kpss_test(df['PMItot'])
KPSS Statistic: 0.2084443853805497
p-value: 0.1
num lags: 9
Critial Values:
   10% : 0.347
   5% : 0.463
   2.5% : 0.574
   1% : 0.739
Result: The series is stationary
In [11]:
kpss_test(df['PMIcomp'])
KPSS Statistic: 0.09140902770821867
p-value: 0.1
num lags: 9
Critial Values:
   10% : 0.347
   5% : 0.463
   2.5% : 0.574
   1% : 0.739
Result: The series is stationary
In [12]:
kpss_test(df['Brent'])
KPSS Statistic: 0.7910625833178427
p-value: 0.01
num lags: 9
Critial Values:
   10% : 0.347
   5% : 0.463
   2.5% : 0.574
   1% : 0.739
Result: The series is not stationary
In [13]:
kpss_test(df['Gold'])
KPSS Statistic: 0.28432725011391896
p-value: 0.1
num lags: 9
Critial Values:
   10% : 0.347
   5% : 0.463
   2.5% : 0.574
   1% : 0.739
Result: The series is stationary
In [14]:
kpss_test(df['SP'])
KPSS Statistic: 1.351825768189659
p-value: 0.01
num lags: 9
Critial Values:
   10% : 0.347
   5% : 0.463
   2.5% : 0.574
   1% : 0.739
Result: The series is not stationary
In [15]:
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['Unemployment']]
ts1['Unemployment'] =ts1['Unemployment'].str.replace('%', '')
ts1=pd.to_numeric(ts1['Unemployment'])
kpss_test(ts1)
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(ts1)
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
KPSS Statistic: 0.3965948161171719
p-value: 0.07862292408742591
num lags: 9
Critial Values:
   10% : 0.347
   5% : 0.463
   2.5% : 0.574
   1% : 0.739
Result: The series is stationary
-2.71677517111288 0.07118014387885119 Ряд не стационарен

Выводы:

Мы проанализировали стационарность и наличие тренда в изучаемых временных рядах. Наши индексы деловой активности явялются стационарными рядами с детерминированным трендом, GDP по тесту ADF также является стационарным детерминированным временным рядом. Особое внимание стоит обратить на уровень безработицы, который имеет нестационарный детерминированный тренд. Индекс S&P 500, фьючерсы на нефть и золото явялются нестационарными рядами со стохастическим временным трендом.

На графике ВВП мы видим характерные черты сезонности, проверим нашу гипотезу. Обнаружение и удаление сезонности.

Есть два способа:

1 Дифференциация: разница в значении во время определенного количества лагов.

2 Декомпозиция: моделирование трендов и сезонности отдельно удаляет их

y(t) = Level + Trend + Seasonality + Noise

In [16]:
from statsmodels.tsa.seasonal import seasonal_decompose
df = pd.read_csv("final.csv", parse_dates=True,sep=';')

df['Date'] = pd.to_datetime(df['Date'])

df.set_index('Date', inplace=True)

df=df[['GDP']]
analysis = df[['GDP']].copy()

decompose_result_additive = seasonal_decompose(analysis, model="ad",period=12)

trend = decompose_result_additive.trend
seasonal = decompose_result_additive.seasonal
residual = decompose_result_additive.resid
decompose_result_additive .plot()
Out[16]:

Как мы видим из графика сезонность есть.В дальнейшем мы будем использовать этот факт для составления прогноза для индексов деловой активности.Проверим наши исследования:

Подбор параметров ARMA- модели.Первые прогнозы:

ARIMA models can be applied only in stationary data.

In [17]:
from statsmodels.tsa.arima_model import ARIMA
import pmdarima as pm
from itertools import product
import pmdarima as pm
from statsmodels.graphics.tsaplots import acf, pacf, plot_acf, plot_pacf
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts=df[['GDP']]
ts['GDP'] = (ts['GDP'] - ts['GDP'].mean()) / ts['GDP'].std() #приведение к нормальному распределению
n_sample = len(ts['GDP'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts['GDP'].iloc[:pivot], ts['GDP'].iloc[pivot:]
In [18]:
arima_model = pm.auto_arima(x_train,start_p=1, start_q=1,
                         test='adf',
                         max_p=3, max_q=3, m=12,
                         start_P=0, seasonal=True,
                         d=None, D=1, trace=True,
                         error_action='ignore',  
                         suppress_warnings=True, 
                         stepwise=True)


arima_model.summary()
plot_pacf(ts).show()
The history saving thread hit an unexpected error (OperationalError('database or disk is full')).History will not be written to the database.
Performing stepwise search to minimize aic
 ARIMA(1,0,1)(0,1,1)[12] intercept   : AIC=inf, Time=0.82 sec
 ARIMA(0,0,0)(0,1,0)[12] intercept   : AIC=411.949, Time=0.01 sec
 ARIMA(1,0,0)(1,1,0)[12] intercept   : AIC=364.371, Time=0.13 sec
 ARIMA(0,0,1)(0,1,1)[12] intercept   : AIC=inf, Time=0.27 sec
 ARIMA(0,0,0)(0,1,0)[12]             : AIC=410.487, Time=0.02 sec
 ARIMA(1,0,0)(0,1,0)[12] intercept   : AIC=389.105, Time=0.04 sec
 ARIMA(1,0,0)(2,1,0)[12] intercept   : AIC=355.260, Time=0.37 sec
 ARIMA(1,0,0)(2,1,1)[12] intercept   : AIC=337.610, Time=0.63 sec
 ARIMA(1,0,0)(1,1,1)[12] intercept   : AIC=inf, Time=0.32 sec
 ARIMA(1,0,0)(2,1,2)[12] intercept   : AIC=inf, Time=1.97 sec
 ARIMA(1,0,0)(1,1,2)[12] intercept   : AIC=inf, Time=1.93 sec
 ARIMA(0,0,0)(2,1,1)[12] intercept   : AIC=358.845, Time=0.53 sec
 ARIMA(2,0,0)(2,1,1)[12] intercept   : AIC=339.519, Time=0.87 sec
 ARIMA(1,0,1)(2,1,1)[12] intercept   : AIC=339.552, Time=0.97 sec
 ARIMA(0,0,1)(2,1,1)[12] intercept   : AIC=341.404, Time=0.86 sec
 ARIMA(2,0,1)(2,1,1)[12] intercept   : AIC=341.142, Time=1.97 sec
 ARIMA(1,0,0)(2,1,1)[12]             : AIC=337.104, Time=0.57 sec
 ARIMA(1,0,0)(1,1,1)[12]             : AIC=335.948, Time=0.23 sec
 ARIMA(1,0,0)(0,1,1)[12]             : AIC=334.029, Time=0.17 sec
 ARIMA(1,0,0)(0,1,0)[12]             : AIC=387.336, Time=0.02 sec
 ARIMA(1,0,0)(0,1,2)[12]             : AIC=335.925, Time=0.41 sec
 ARIMA(1,0,0)(1,1,0)[12]             : AIC=362.654, Time=0.06 sec
 ARIMA(1,0,0)(1,1,2)[12]             : AIC=inf, Time=1.57 sec
 ARIMA(0,0,0)(0,1,1)[12]             : AIC=358.038, Time=0.11 sec
 ARIMA(2,0,0)(0,1,1)[12]             : AIC=335.944, Time=0.25 sec
 ARIMA(1,0,1)(0,1,1)[12]             : AIC=335.971, Time=0.26 sec
 ARIMA(0,0,1)(0,1,1)[12]             : AIC=338.511, Time=0.18 sec
 ARIMA(2,0,1)(0,1,1)[12]             : AIC=338.029, Time=0.25 sec
 ARIMA(1,0,0)(0,1,1)[12] intercept   : AIC=inf, Time=0.30 sec

Best model:  ARIMA(1,0,0)(0,1,1)[12]          
Total fit time: 16.118 seconds
In [19]:
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['Gold']]
ts1['Gold'] = (ts1['Gold'] - ts1['Gold'].mean()) / ts1['Gold'].std() #приведение к нормальному распределению
n_sample = len(ts1['Gold'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['Gold'].iloc[:pivot], ts1['Gold'].iloc[pivot:]
arima_model = pm.auto_arima(x_train, 
                              start_p=1,d =0,start_q=0,
                              test="adf", supress_warnings = True,
                              trace=True)
arima_model.summary()
plot_pacf(ts1).show()
Performing stepwise search to minimize aic
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=123.644, Time=0.04 sec
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=420.585, Time=0.01 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=283.501, Time=0.06 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=420.050, Time=0.00 sec
 ARIMA(2,0,0)(0,0,0)[0] intercept   : AIC=124.516, Time=0.05 sec
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=124.278, Time=0.09 sec
 ARIMA(2,0,1)(0,0,0)[0] intercept   : AIC=127.644, Time=0.10 sec
 ARIMA(1,0,0)(0,0,0)[0]             : AIC=121.645, Time=0.05 sec
 ARIMA(2,0,0)(0,0,0)[0]             : AIC=122.518, Time=0.06 sec
 ARIMA(1,0,1)(0,0,0)[0]             : AIC=122.279, Time=0.04 sec
 ARIMA(0,0,1)(0,0,0)[0]             : AIC=282.567, Time=0.03 sec
 ARIMA(2,0,1)(0,0,0)[0]             : AIC=123.567, Time=0.10 sec

Best model:  ARIMA(1,0,0)(0,0,0)[0]          
Total fit time: 0.652 seconds
In [20]:
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['PMItot']]
ts1['PMItot'] = (ts1['PMItot'] - ts1['PMItot'].mean()) / ts1['PMItot'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMItot'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMItot'].iloc[:pivot], ts1['PMItot'].iloc[pivot:]
arima_model = pm.auto_arima(x_train, 
                              start_p=1,d =0,start_q=0,
                              test="adf", supress_warnings = True,
                              trace=True)
arima_model.summary()
plot_pacf(ts1).show()
Performing stepwise search to minimize aic
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=188.288, Time=0.03 sec
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=350.898, Time=0.02 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=248.906, Time=0.05 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=358.367, Time=0.01 sec
 ARIMA(2,0,0)(0,0,0)[0] intercept   : AIC=188.598, Time=0.06 sec
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=188.645, Time=0.06 sec
 ARIMA(2,0,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.22 sec
 ARIMA(1,0,0)(0,0,0)[0]             : AIC=187.389, Time=0.01 sec
 ARIMA(2,0,0)(0,0,0)[0]             : AIC=187.886, Time=0.03 sec
 ARIMA(1,0,1)(0,0,0)[0]             : AIC=187.908, Time=0.03 sec
 ARIMA(0,0,1)(0,0,0)[0]             : AIC=253.727, Time=0.02 sec
 ARIMA(2,0,1)(0,0,0)[0]             : AIC=189.032, Time=0.06 sec

Best model:  ARIMA(1,0,0)(0,0,0)[0]          
Total fit time: 0.624 seconds
In [21]:
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['PMIcomp']]
ts1['PMIcomp'] = (ts1['PMIcomp'] - ts1['PMIcomp'].mean()) / ts1['PMIcomp'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMIcomp'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMIcomp'].iloc[:pivot], ts1['PMIcomp'].iloc[pivot:]
arima_model = pm.auto_arima(x_train, 
                              start_p=1,d =0,start_q=0,
                              test="adf", supress_warnings = True,
                              trace=True)
arima_model.summary()
plot_pacf(ts1).show()
Performing stepwise search to minimize aic
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=106.151, Time=0.03 sec
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=204.512, Time=0.02 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=148.596, Time=0.03 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=203.761, Time=0.02 sec
 ARIMA(2,0,0)(0,0,0)[0] intercept   : AIC=105.420, Time=0.05 sec
 ARIMA(3,0,0)(0,0,0)[0] intercept   : AIC=107.283, Time=0.11 sec
 ARIMA(2,0,1)(0,0,0)[0] intercept   : AIC=107.096, Time=0.09 sec
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=105.623, Time=0.08 sec
 ARIMA(3,0,1)(0,0,0)[0] intercept   : AIC=108.822, Time=0.20 sec
 ARIMA(2,0,0)(0,0,0)[0]             : AIC=103.429, Time=0.03 sec
 ARIMA(1,0,0)(0,0,0)[0]             : AIC=104.203, Time=0.02 sec
 ARIMA(3,0,0)(0,0,0)[0]             : AIC=105.299, Time=0.06 sec
 ARIMA(2,0,1)(0,0,0)[0]             : AIC=105.114, Time=0.06 sec
 ARIMA(1,0,1)(0,0,0)[0]             : AIC=103.629, Time=0.03 sec
 ARIMA(3,0,1)(0,0,0)[0]             : AIC=106.835, Time=0.12 sec

Best model:  ARIMA(2,0,0)(0,0,0)[0]          
Total fit time: 0.955 seconds
In [22]:
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['USD-EUR']]
ts1['USD-EUR'] = (ts1['USD-EUR'] - ts1['USD-EUR'].mean()) / ts1['USD-EUR'].std() #приведение к нормальному распределению
n_sample = len(ts1['USD-EUR'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['USD-EUR'].iloc[:pivot], ts1['USD-EUR'].iloc[pivot:]
arima_model = pm.auto_arima(x_train, 
                              start_p=1,d =0,start_q=0,
                              test="adf", supress_warnings = True,
                              trace=True)
arima_model.summary()
plot_pacf(ts1).show()
Performing stepwise search to minimize aic
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=inf, Time=0.03 sec
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=356.967, Time=0.02 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=196.537, Time=0.04 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=394.848, Time=0.01 sec
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=-155.775, Time=0.11 sec
 ARIMA(2,0,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.21 sec
 ARIMA(1,0,2)(0,0,0)[0] intercept   : AIC=-153.778, Time=0.18 sec
 ARIMA(0,0,2)(0,0,0)[0] intercept   : AIC=86.210, Time=0.08 sec
 ARIMA(2,0,0)(0,0,0)[0] intercept   : AIC=inf, Time=0.15 sec
 ARIMA(2,0,2)(0,0,0)[0] intercept   : AIC=-152.119, Time=0.36 sec
 ARIMA(1,0,1)(0,0,0)[0]             : AIC=-157.661, Time=0.09 sec
 ARIMA(0,0,1)(0,0,0)[0]             : AIC=229.348, Time=0.04 sec
 ARIMA(1,0,0)(0,0,0)[0]             : AIC=inf, Time=0.03 sec
 ARIMA(2,0,1)(0,0,0)[0]             : AIC=-155.545, Time=0.07 sec
 ARIMA(1,0,2)(0,0,0)[0]             : AIC=-155.664, Time=0.14 sec
 ARIMA(0,0,2)(0,0,0)[0]             : AIC=111.458, Time=0.12 sec
 ARIMA(2,0,0)(0,0,0)[0]             : AIC=inf, Time=0.08 sec
 ARIMA(2,0,2)(0,0,0)[0]             : AIC=-154.048, Time=0.23 sec

Best model:  ARIMA(1,0,1)(0,0,0)[0]          
Total fit time: 1.999 seconds
In [23]:
#приведем данные к стационарному виду
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
df['Unemployment'] = df['Unemployment'].str.replace('%', '')
df['Unemployment']=pd.to_numeric(df['Unemployment'])
df['Unemployment_dif']=(df['Unemployment']-df['Unemployment'].shift(1))
df[['Unemployment_dif']].dropna()
df= df[df['Unemployment_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Unemployment_dif'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-11.432120662348963 6.478415547715242e-21 Ряд стационарен
In [24]:
#приведем данные к стационарному виду
df['Brent_dif']=(df['Brent']-df['Brent'].shift(1)).dropna()
df[['Brent_dif']].dropna()
df= df[df['Brent_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Brent_dif'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-6.3209632603841 3.06699569955005e-08 Ряд стационарен
In [25]:
df['Gold_dif']=(df['Gold']-df['Gold'].shift(1)).dropna()
df[['Gold_dif']].dropna()
df= df[df['Gold_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Gold_dif'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-11.35837278788237 9.608487353219511e-21 Ряд стационарен
In [26]:
df['SP_dif']=(df['SP']-df['SP'].shift(1)).dropna()
df[['SP_dif']].dropna()
df= df[df['SP_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['SP_dif'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-5.647512295305509 1.003450715228814e-06 Ряд стационарен
In [27]:
df['USD-EUR_dif']=(df['USD-EUR']-df['USD-EUR'].shift(1)).dropna()
df[['USD-EUR_dif']].dropna()
df= df[df['USD-EUR_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['USD-EUR_dif'])
stat, pvalue
if pvalue < alpha:
    print(stat, pvalue,'Ряд стационарен')
else:
    print(stat, pvalue,'Ряд не стационарен')
-15.018877771046922 1.0293226979823363e-27 Ряд стационарен

Первые прогнозы

In [28]:
import statsmodels.api as sm
import pmdarima as pm
from itertools import product
import pmdarima as pm
#df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['PMItot']]
ts1['PMItot'] = (ts1['PMItot'] - ts1['PMItot'].mean()) / ts1['PMItot'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMItot'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMItot'].iloc[:pivot], ts1['PMItot'].iloc[pivot:]
first_val = x_train.iloc[0]
x_train_diff = x_train.diff().dropna()
model = sm.tsa.arima.ARIMA(x_train_diff.values, order=(1,0,0)).fit()
x_pred_diff = model.forecast(len(x_test) - 1)
x_pred_diff = np.array([first_val] + x_pred_diff.tolist())
x_pred = pd.Series(np.cumsum(x_pred_diff))
plt.plot(np.arange(len(x_pred)) + len(x_train),x_test)
plt.plot(np.arange(len(x_pred)) + len(x_train), x_pred)
plt.plot(x_train)
plt.plot(x_train)

plt.legend(['y_test', 'y_test_pred', 'y_train','y_train_pred'])
Out[28]:
<matplotlib.legend.Legend at 0x1cb37666a48>
In [29]:
ts1=df[['PMIcomp']]
ts1['PMIcomp'] = (ts1['PMIcomp'] - ts1['PMIcomp'].mean()) / ts1['PMIcomp'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMIcomp'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMIcomp'].iloc[:pivot], ts1['PMIcomp'].iloc[pivot:]
first_val = x_train.iloc[0]
x_train_diff = x_train.diff().dropna()
model = sm.tsa.arima.ARIMA(x_train_diff.values, order=(2,0,0)).fit()
x_pred_diff = model.forecast(len(x_test)-1 )
x_pred_diff = np.array([first_val] + x_pred_diff.tolist())
x_pred = pd.Series(np.cumsum(x_pred_diff))
plt.plot(np.arange(len(x_pred)) + len(x_train),x_test)
plt.plot(np.arange(len(x_pred)) + len(x_train), x_pred)
plt.plot(x_train)
plt.plot(x_train)

plt.legend(['y_test', 'y_test_pred', 'y_train','y_train_pred'])
Out[29]:
<matplotlib.legend.Legend at 0x1cb346a3f08>

Как мы видим, несмотря на то, что модели, подобранные алгоритмом, максимизируют точность - актуальный прогноз построить не получилось, поэтому расмотрим подробнее зависимые переменые:

Множественная регрессия временных рядов:

Общая регрессионная модель временных рядов с несколькими объясняющими переменными является расширением авторегрессионной модели с распределенными лагами на случай включения в число регрессоров запаздываний нескольких переменных. В общем случае регрессионная модель временных рядов расширяется на случай k дополнительных регрессов и включает q1:

  • запаздывания первой объясняющей переменной, q1
  • запаздывания второй объясняющей переменной и так далее.

Однако большинство наших временных рядов являются нестационнарными временными рядами.И прежде чем составлять какой-либо моножественный прогноз мы должны выяснить

1) Обладают ли ряды детерминированными или стохастическими трендами

2) Устранение нестационарности: "на практике для стационарности многомерного ряда достаточно стационарности составляющих его одномерных рядов."

3) Провести анализ на наличие структурных разрывов

4) Провести тест на причинность по Гренджеру

4) Проверить коинтегрированность индексов в случае отсутствия ложной регрессии между ними (тест ADF-GL/ тест Йенсена)

5) Построить VAR() модель : оценить параметры; Построить многомерный прогноз;

Cross Correlation Function (CCF) and ADF-GL test + Johansen test for co-integration

In [30]:
# CFF-test and Granger-Causality test
cff=sm.tsa.stattools.ccf(df['SP_dif'],df['PMItot'], adjusted=False)
cff[:7]
Out[30]:
array([ 0.07488947,  0.0553094 ,  0.02316356, -0.01985263,  0.04202412,
        0.04428844,  0.04700216])
In [31]:
cff=sm.tsa.stattools.ccf(df['Gold_dif'],df['PMItot'], adjusted=False)
cff[:7]
Out[31]:
array([0.22606771, 0.16937807, 0.17561882, 0.12993433, 0.17316124,
       0.14612628, 0.10260552])
In [32]:
cff=sm.tsa.stattools.ccf(df['Brent_dif'],df['PMItot'], adjusted=False)
cff[:7]
Out[32]:
array([ 0.04014853,  0.0378068 ,  0.02146401, -0.0099423 ,  0.01327634,
       -0.0325728 , -0.02874329])
In [33]:
cff=sm.tsa.stattools.ccf(df['GDP'],df['PMItot'], adjusted=False)
cff[:7]
Out[33]:
array([0.22139378, 0.21931753, 0.24745512, 0.27972358, 0.30801513,
       0.19946442, 0.0719124 ])
In [34]:
cff=sm.tsa.stattools.ccf(df['USD-EUR'],df['PMItot'], adjusted=False)
cff[:7]
Out[34]:
array([0.21253023, 0.19421706, 0.17772199, 0.16501161, 0.14197805,
       0.12735322, 0.10884041])
In [35]:
cff=sm.tsa.stattools.ccf(df['Unemployment'],df['PMItot'], adjusted=False)
cff[:7]
Out[35]:
array([-0.03933208, -0.1294871 , -0.15172892, -0.157741  , -0.15299171,
       -0.17251836, -0.18594392])
In [36]:
cff=sm.tsa.stattools.ccf(df['SP_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
Out[36]:
array([-0.03426587, -0.03967198, -0.10715686, -0.07437971, -0.0222114 ,
        0.000311  ,  0.02885604])
In [37]:
cff=sm.tsa.stattools.ccf(df['Gold_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
Out[37]:
array([0.01830157, 0.01874825, 0.04640013, 0.03670187, 0.06211483,
       0.04169818, 0.08650913])
In [38]:
cff=sm.tsa.stattools.ccf(df['Brent_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
Out[38]:
array([-0.02361618, -0.03545415,  0.03700597,  0.00578842, -0.01281549,
       -0.05511056, -0.0946557 ])
In [39]:
cff=sm.tsa.stattools.ccf(df['GDP'],df['PMIcomp'], adjusted=False)
cff[:7]
Out[39]:
array([-0.05335246, -0.05874772, -0.02140082,  0.00806623,  0.07130947,
        0.04800688,  0.01377704])
In [40]:
cff=sm.tsa.stattools.ccf(df['Unemployment_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
Out[40]:
array([ 0.01963671, -0.02212346, -0.03394327, -0.02038006, -0.02043086,
       -0.01232782, -0.01546808])
In [41]:
cff=sm.tsa.stattools.ccf(df['USD-EUR'],df['PMIcomp'], adjusted=False)
cff[:7]
Out[41]:
array([-0.0642212 , -0.07644561, -0.08435376, -0.08392729, -0.08074017,
       -0.08666067, -0.09400758])

Выводы:

Тест на причинность по Гренджеру T he p-value of each test is less than 5%, so it can be said that the any index is useful for forecasting Index of intepreners' activity .

In [42]:
from statsmodels.tsa.stattools import grangercausalitytests
grangercausalitytests(df[['PMItot', 'USD-EUR']], maxlag=[2])
Granger Causality
number of lags (no zero) 2
ssr based F test:         F=1.2414  , p=0.2913  , df_denom=194, df_num=2
ssr based chi2 test:   chi2=2.5467  , p=0.2799  , df=2
likelihood ratio test: chi2=2.5306  , p=0.2822  , df=2
parameter F test:         F=1.2414  , p=0.2913  , df_denom=194, df_num=2
Out[42]:
{2: ({'ssr_ftest': (1.2413689663134728, 0.29127343627097546, 194.0, 2),
   'ssr_chi2test': (2.5467260236740317, 0.27988876801294515, 2),
   'lrtest': (2.5305677224930605, 0.28215919071205225, 2),
   'params_ftest': (1.2413689663133483, 0.2912734362710129, 194.0, 2.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb1f259488>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb1db17f48>,
   array([[0., 0., 1., 0., 0.],
          [0., 0., 0., 1., 0.]])])}
In [43]:
grangercausalitytests(df[['PMItot', 'GDP']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=2.0602  , p=0.1528  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=2.0916  , p=0.1481  , df=1
likelihood ratio test: chi2=2.0807  , p=0.1492  , df=1
parameter F test:         F=2.0602  , p=0.1528  , df_denom=197, df_num=1
Out[43]:
{1: ({'ssr_ftest': (2.060220088874574, 0.15277412703233398, 197.0, 1),
   'ssr_chi2test': (2.0915939988574355, 0.14811146248394463, 1),
   'lrtest': (2.0807327440405743, 0.14916855428311235, 1),
   'params_ftest': (2.0602200888746007, 0.15277412703233398, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2f771408>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb1ba7c0c8>,
   array([[0., 1., 0.]])])}
In [44]:
grangercausalitytests(df[['PMItot','Gold_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.1205  , p=0.7289  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.1223  , p=0.7265  , df=1
likelihood ratio test: chi2=0.1223  , p=0.7266  , df=1
parameter F test:         F=0.1205  , p=0.7289  , df_denom=197, df_num=1
Out[44]:
{1: ({'ssr_ftest': (0.12049868176489817, 0.7288648059690301, 197.0, 1),
   'ssr_chi2test': (0.12233368707096262, 0.7265170722588772, 1),
   'lrtest': (0.12229628849308938, 0.7265572018703386, 1),
   'params_ftest': (0.12049868176501637, 0.7288648059689118, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2e767408>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2e767c48>,
   array([[0., 1., 0.]])])}
In [45]:
grangercausalitytests(df[['PMItot','Brent_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.5589  , p=0.4556  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.5674  , p=0.4513  , df=1
likelihood ratio test: chi2=0.5666  , p=0.4516  , df=1
parameter F test:         F=0.5589  , p=0.4556  , df_denom=197, df_num=1
Out[45]:
{1: ({'ssr_ftest': (0.5589313469762955, 0.4555826145051268, 197.0, 1),
   'ssr_chi2test': (0.5674429918541071, 0.45127679609768534, 1),
   'lrtest': (0.5666395323478355, 0.4515973752845611, 1),
   'params_ftest': (0.5589313469763506, 0.4555826145051003, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb32ecf3c8>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb32ecf808>,
   array([[0., 1., 0.]])])}
In [46]:
grangercausalitytests(df[['PMItot','SP_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.1764  , p=0.6749  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.1791  , p=0.6722  , df=1
likelihood ratio test: chi2=0.1790  , p=0.6722  , df=1
parameter F test:         F=0.1764  , p=0.6749  , df_denom=197, df_num=1
Out[46]:
{1: ({'ssr_ftest': (0.17640140527268897, 0.6749419462642945, 197.0, 1),
   'ssr_chi2test': (0.17908772108902435, 0.6721584135533447, 1),
   'lrtest': (0.17900758789210158, 0.672227494537481, 1),
   'params_ftest': (0.17640140527279857, 0.6749419462641995, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb32ed1488>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb32ed1a08>,
   array([[0., 1., 0.]])])}
In [47]:
grangercausalitytests(df[['PMIcomp','SP_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.0304  , p=0.8618  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.0308  , p=0.8606  , df=1
likelihood ratio test: chi2=0.0308  , p=0.8606  , df=1
parameter F test:         F=0.0304  , p=0.8618  , df_denom=197, df_num=1
Out[47]:
{1: ({'ssr_ftest': (0.030368369812070573, 0.8618359779890338, 197.0, 1),
   'ssr_chi2test': (0.030830832296518347, 0.8606183317622438, 1),
   'lrtest': (0.03082845619019281, 0.8606236479107884, 1),
   'params_ftest': (0.030368369812047116, 0.861835977989083, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb32ecf6c8>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2e767cc8>,
   array([[0., 1., 0.]])])}
In [48]:
grangercausalitytests(df[['PMIcomp','Brent_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.1360  , p=0.7127  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.1381  , p=0.7102  , df=1
likelihood ratio test: chi2=0.1380  , p=0.7102  , df=1
parameter F test:         F=0.1360  , p=0.7127  , df_denom=197, df_num=1
Out[48]:
{1: ({'ssr_ftest': (0.1360199381735314, 0.7126662927530032, 197.0, 1),
   'ssr_chi2test': (0.1380913077903872, 0.7101859320565205, 1),
   'lrtest': (0.13804365670000607, 0.710233680171941, 1),
   'params_ftest': (0.13601993817351837, 0.7126662927530032, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb3766c308>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb3766c108>,
   array([[0., 1., 0.]])])}
In [49]:
grangercausalitytests(df[['PMIcomp','GDP']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.6846  , p=0.4090  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.6950  , p=0.4045  , df=1
likelihood ratio test: chi2=0.6938  , p=0.4049  , df=1
parameter F test:         F=0.6846  , p=0.4090  , df_denom=197, df_num=1
Out[49]:
{1: ({'ssr_ftest': (0.6846079499038269, 0.4090056356853612, 197.0, 1),
   'ssr_chi2test': (0.6950334516790121, 0.4044575757294545, 1),
   'lrtest': (0.6938285635832244, 0.40486518918896475, 1),
   'params_ftest': (0.684607949903744, 0.4090056356853913, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2f7760c8>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2f776808>,
   array([[0., 1., 0.]])])}
In [50]:
grangercausalitytests(df[['PMIcomp','Unemployment_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.2987  , p=0.5853  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.3033  , p=0.5818  , df=1
likelihood ratio test: chi2=0.3031  , p=0.5820  , df=1
parameter F test:         F=0.2987  , p=0.5853  , df_denom=197, df_num=1
Out[50]:
{1: ({'ssr_ftest': (0.29874228332809166, 0.585290496084301, 197.0, 1),
   'ssr_chi2test': (0.30329165820110826, 0.5818261676859093, 1),
   'lrtest': (0.30306192584998826, 0.5819692052982772, 1),
   'params_ftest': (0.29874228332806774, 0.5852904960843148, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb346b32c8>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb346b3488>,
   array([[0., 1., 0.]])])}
In [51]:
grangercausalitytests(df[['PMIcomp','Gold_dif']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.1009  , p=0.7511  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=0.1025  , p=0.7489  , df=1
likelihood ratio test: chi2=0.1024  , p=0.7489  , df=1
parameter F test:         F=0.1009  , p=0.7511  , df_denom=197, df_num=1
Out[51]:
{1: ({'ssr_ftest': (0.10091908380940975, 0.7510659393984734, 197.0, 1),
   'ssr_chi2test': (0.10245592264914696, 0.7489021017980605, 1),
   'lrtest': (0.10242968856800871, 0.748933168132953, 1),
   'params_ftest': (0.10091908380934857, 0.7510659393985517, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2a08b788>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb2a08b748>,
   array([[0., 1., 0.]])])}
In [52]:
grangercausalitytests(df[['PMIcomp','PMItot']], maxlag=[1])
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=4.2859  , p=0.0397  , df_denom=197, df_num=1
ssr based chi2 test:   chi2=4.3511  , p=0.0370  , df=1
likelihood ratio test: chi2=4.3045  , p=0.0380  , df=1
parameter F test:         F=4.2859  , p=0.0397  , df_denom=197, df_num=1
Out[52]:
{1: ({'ssr_ftest': (4.285875254184534, 0.03973437006593587, 197.0, 1),
   'ssr_chi2test': (4.35114239003506, 0.03698412989998138, 1),
   'lrtest': (4.30448676129447, 0.038011963404940785, 1),
   'params_ftest': (4.285875254184432, 0.03973437006593837, 197.0, 1.0)},
  [<statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb1f23fec8>,
   <statsmodels.regression.linear_model.RegressionResultsWrapper at 0x1cb1f23fa88>,
   array([[0., 1., 0.]])])}

Выводы

In [72]:
from statsmodels.tsa.vector_ar.vecm import coint_johansen
from statsmodels.tsa.vector_ar.vecm import coint_johansen
 # dataframe of n series for cointegration analysis


import numpy as np
from statsmodels.tsa.tsatools import lagmat
 
 
def johansen(ts, lags):
    """
    Calculate the Johansen Test for the given time series
    """
    # Make sure we are working with arrays, convert if necessary
    ts = np.asarray(ts)
 
    # nobs is the number of observations
    # m is the number of series
    nobs, m = ts.shape
 
    # Obtain the cointegrating vectors and corresponding eigenvalues
    eigenvectors, eigenvalues = maximum_likelihood_estimation(ts, lags)
 
    # Build table of critical values for the hypothesis test
    critical_values_string = """2.71   3.84    6.63
             13.43   15.50   19.94
             27.07   29.80   35.46
             44.49   47.86   54.68
             65.82   69.82   77.82
             91.11   95.75   104.96
             120.37  125.61  135.97
             153.63  159.53  171.09
             190.88  197.37  210.06
             232.11  239.25  253.24
             277.38  285.14  300.29
             326.53  334.98  351.25"""
    select_critical_values = np.array(
            critical_values_string.split(),
            float).reshape(-1, 3)
    critical_values = select_critical_values[:, 1]
 
    # Calculate numbers of cointegrating relations for which
    # the null hypothesis is rejected
    rejected_r_values = []
    for r in range(m):
        if h_test(eigenvalues, r, nobs, lags, critical_values):
            rejected_r_values.append(r)
 
    return eigenvectors, rejected_r_values

 
def h_test(eigenvalues, r, nobs, lags, critical_values):
    """
    Helper to execute the hypothesis test
    """
    # Calculate statistic
    t = nobs - lags - 1
    m = len(eigenvalues)
    statistic = -t * np.sum(np.log(np.ones(m) - eigenvalues)[r:])
 
    # Get the critical value
    critical_value = critical_values[m - r - 1]
 
    # Finally, check the hypothesis
    if statistic > critical_value:
        return True
    else:
        return False
 
 
def maximum_likelihood_estimation(ts, lags):
    """
    Obtain the cointegrating vectors and corresponding eigenvalues
    """
    # Make sure we are working with array, convert if necessary
    ts = np.asarray(ts)
 
    # Calculate the differences of ts
    ts_diff = np.diff(ts, axis=0)
 
    # Calculate lags of ts_diff.
    ts_diff_lags = lagmat(ts_diff, lags, trim='both')
 
    # First lag of ts
    ts_lag = lagmat(ts, 1, trim='both')
 
    # Trim ts_diff and ts_lag
    ts_diff = ts_diff[lags:]
    ts_lag = ts_lag[lags:]
 
    # Include intercept in the regressions
    ones = np.ones((ts_diff_lags.shape[0], 1))
    ts_diff_lags = np.append(ts_diff_lags, ones, axis=1)
 
    # Calculate the residuals of the regressions of diffs and lags
    # into ts_diff_lags
    inverse = np.linalg.pinv(ts_diff_lags)
    u = ts_diff - np.dot(ts_diff_lags, np.dot(inverse, ts_diff))
    v = ts_lag - np.dot(ts_diff_lags, np.dot(inverse, ts_lag))
 
    # Covariance matrices of the calculated residuals
    t = ts_diff_lags.shape[0]
    Svv = np.dot(v.T, v) / t
    Suu = np.dot(u.T, u) / t
    Suv = np.dot(u.T, v) / t
    Svu = Suv.T
 
    # ToDo: check for singular matrices and exit
    Svv_inv = np.linalg.inv(Svv)
    Suu_inv = np.linalg.inv(Suu)
 
    # Calculate eigenvalues and eigenvectors of the product of covariances
    cov_prod = np.dot(Svv_inv, np.dot(Svu, np.dot(Suu_inv, Suv)))
    eigenvalues, eigenvectors = np.linalg.eig(cov_prod)
 
    # Use Cholesky decomposition on eigenvectors
    evec_Svv_evec = np.dot(eigenvectors.T, np.dot(Svv, eigenvectors))
    cholesky_factor = np.linalg.cholesky(evec_Svv_evec)
    eigenvectors = np.dot(eigenvectors, np.linalg.inv(cholesky_factor.T))
 
    # Order the eigenvalues and eigenvectors
    indices_ordered = np.argsort(eigenvalues)
    indices_ordered = np.flipud(indices_ordered)
 
    # Return the calculated values
    return eigenvalues[indices_ordered], eigenvectors[:, indices_ordered]

johansen(df[['Brent','Gold','PMItot','PMIcomp','USD-EUR','SP','GDP']],3)
Out[72]:
(array([0.31079533, 0.15707143, 0.12569789, 0.09625334, 0.07158033,
        0.02363577, 0.00067435]), [])
In [ ]: